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ABSTRACT 


The objective is to obtain and compare interpolation 
methods and contour plot techniques, based on data from 
Sites whose locations are irregular or scattered. Data 
poem the 1980 Fugro report (ieee. 1] on MX valley soil 
samples serves as a test bed. 

Two groups of data, seismic p-wave velocities and sur- 
face soil depth were studied using six different interpola- 
tion methods and three different contour plot techniques. 
Comparison of the interpolation methods and of the contour 


plot techniques are made. 
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I. INTRODUCTION 


Detailed investigations of seismic p-wave velocities and 
depth of the surface soil layer were performed in the Ralston 
Valley, Nevada, by Fugro, ines Ret . el Qi enews a oA em, 
Engineer Waterways Experiment Station (W.E.S.) as part of the 
preliminary work on the MX missile system deployment. Six- 
teen locations were chosen in the Ralston Valley by W.E.S. 
for their study organized into four sets of sites (coded RA, 
RB,RC,RD) located in each of the four predominant surficial 
soll types (coded 51,5Y,U,4U), and an additional site marked 
RU2 identified in the Valley. 

The data source for identified sites [Ref. 2] gives the 
(X,Y) coordinates (one mile to the inch) of the boring 
locations relative to the map's lower left grid mark, the 
seismic p-wave velocities (MPS - meters per second) of the 
Surface layer, and the depth of that layer (meters) for the 
2B line; see Appendix A. The site identifiers and locations 
for the Ralston Valley sampling plan can be found in Appendix 
B. The immediate concern is to develop contours describing 
the depth of the top layer, and of the compression wave speed 
in the top layer, in the entire Valley. 

Based on these limited data, the mathematical problem is 
that of constructing a smooth bivariate function F(X,Y) which 


takes on specified values F(X ,Y,) = PL omy w= eh ateg Ne 


2 
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Franke [ Ref. Bp), 5 AOL has made a rather extensive general 

study of the modern interpolation methods. Most of the methods 
require the user to specify a parameter (or several), but 

these parameters can be determined by computational experi- 
ment [Ref. 4]. The interpolation methods selected for 
Seeeittc inelusion im the current work are a subset of these 
methods, see Table 1. Furthermore, contour plotting, specific~ 
ally non-IMSL subroutine CONTUR, subroutine CONISD and sub- 
routine PLT3Dl, were exercised for comparison. Basic FORTRAN 
programs for the first five interpolation methods were obtained 
from Professor R. Franke [ Ref. By475, 60! « 

The various interpolation methods were uSed to supply 
function values at a gridwork of reference points. In order 
to compare the methods, the mean and standard deviation of 
the values were produced at each grid point. Also it is 
important to identify which methods generates sharply different 
interpolation values compared to the other methods. For this 
purpose, the same FORTRAN program generates the absolute dif- 
ferences from the mean values for each method. 

FORTRAN programs for all of the methods can be obtained 
by writing the advisor. 

In the next chapter we will describe a number of inter- 
polation methods. Chapter 3 describes employment of the 
Ralston Valley data as a way of comparing these methods. 
Chapter 4 discusses contour plot techniques, and the con- 
clusions from using these procedures are given in Chapter 5. 


Graphical results are presented in a series of appendices. 


os 








a ia. —— 
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3) 


4) 


5) 


6) 


TABLE lL 
INTERPOLATION METHODS 


Inverse Distance Weighted Methods 

A) Shepard's Method 

B) Modified Shepard's Method 

C) Modified Linear Shepard's Method 

D) Modified Shepard's Method Boolean Sum Plane 
E) Modified McLain Method 

F) Quadratic Shepard's Method 
*G) Modified Quadratic Shepard's Method 

H) McLain Method 

Franke's Methods (Rectangle based blending methods) 
A) Franke's Method (Mode One) 

B) Franke's Method (Mode Three) 
“@) Fxranke's Method (Thin plate local functions) 
Triangle Based Blending Methods 

A) Nielson-Franke Linear Triangle Method 
*B) Nielson-Franke Quadratic Triangle Method 
Finite Element Based Methods 

A) Akima's Method 

B) Akima's Method - Modification Two 

C) Akima's Method - Modification One 

D) Akima's Method - Modification Three 
*E) Nielson's Minimum Norm Network 

F) Lawson's Method 

Foley's Method 

A) Generalized Newton Interpolant 

B) TF Delta Sum Berstein Interpolant 

GC)” Teerated Delta Sums: TF Delta Sum Bicubic Spline 
*D) Iterated Delta Sums: A Shepard's Method Delta Sum 

Buc. eS: 

Nodal Basis Function Type Methods 

A) Rotated Gaussian 

*B) Hardy's Multiquadric 

C) Hardy's Reciprocal Multiquadric 

BD) Dueton's Radial Cubic Method 

E) Dwehom'’s Thin Plate Function 

F) Rotated B-Spline 


* Methods chosen for comparison study. 
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It. INTERPOLATION METHODS 


The problem being addressed here is that of constructing 
a smooth bivariate function, F(X,Y), which takes on specified 
values P(X) Vy) = Pus k=1,...,N based on point data from 
Scattered, lrregular locations. The bivariate function, 

Pim ,Y) 1S smooth, (i.e., it has at least continuous first 
partial derivatives), and the points (XY) are 'scattered', 
that 1s, not necessarily conforming to some regular pattern. 

Some interpolation methods, called 'global', are sensitive 
to all the data points: the addition or deletion of a data 
point, or a change of one of the coordinates of a data point, 
W1ll propogate throughout the entire domain of definition. 
Other interpolation methods are called ‘local', that is, insen- 
Sitive to changes in the data points which are sufficiently 
distant from the location at hand. 

Franke, R. [| Ref. 3,4,5,6 ] has made an extensive general 
study of new interpolation methods. He produced a classifi- 
cation of the interpolation methods [Ref. 4], as shown in 
Table 1. Under his six main types, there are a total of 29 
Gifferent methods. All the methods appearing in the same 
group are either local or global and use the same basic idea; 
several are modifications of the previously described ones. 

Six methods, one from each of the six main types, are 


described below. 


Lo 





A. MODIFIED QUADRATIC SHEPARD'S METHOD 


The modified quadratic Shepard's method was derived 
from the inverse distance weighted method. All methods of 
this type which are considered may be viewed as generaliza- 


tions of Shepard's method. The basic Shepard's method is 


N N 
Peep) = 9) 2 We (X,Y).0,.(X,Y)/ © W,. (X,Y) 
k=1 * i k=1 * 


where W,. (X,Y) =d Typically uw = 2, although other values 


k 


Pay. y Z as 2 
may be used Here qd, = ((X Ky) + (Y Y.) ys 


In the modified quadratic Shepard's method the weights 


for obtaining the nodal functions (quadratics) are taken as 


2 
( = ) 
W, (X,Y) = Pe oe and R = 1/2 \[N_/N:D 
s Gee 4.) 2 : 


where D is thediameter of the point set, and a nominal value 


for Nge determined by computational experiment | Ref. 4], 1s 


- = 18. The diameter of the point set is defined 
2 Z Z 
D = MAX (X. - X.) + (Y. - Y=) 
Tee: ele J Ee J 


The nodal functions Q, (X,Y) are defined as 
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a _— 2 


(X,Y) =F + X-X = > 
Q. ie aw * + ae ae ~ af (X ae 


le2 i 

a ES 2 
+ a (X-X (y= YV +r YY 

oe = 7 a a ~ 


and for the coefficients one must solve the least squares 


problem, 
N 
MIN UN F + = _- 
fa, fF, eee ee 
Bg mon’ diaghk 


2 

> Ag (Xi -X,) + ee ee 
2 2 

- ag (5 -¥,) - | 


Complete details are given in | Ref. See 


me LOCAL THIN PLATE SPLINES METHOD 
This is known as Franke's method [ Ref. 4 |. The idea is 
to represent the interpolation function as 
re) Wet xX,Y)°O. . (X,Y) 
i,j 4) 1) 
where the oe) are local weight functions | Ref. 6 |. Here 
we choose the W_ so that ec 1, and the local approximations 


1L- 
1) 
os 5 interpolate at points where the corresponding weight 


function is non-zero. The weight functions are products of 


plecewise hermite cubics [ Ref. 6 |. 


iy 





The local approximations for this method are taken to be 


the thin plate splines and have the form 


2 
Oey = 5 Ad “Pee a + at b*X + GY 
td ker * k 5 
2 : 2 
where dv = ((X-X,) + (¥-¥,)~) 


Mma I is the set of indices k for which (XY). FL) 1S a point 


to be interpolated by Q15°%,¥), and is defined by 
I= { x >: QO 1S to take the value FL at 57%.) } 


The coefficients A., and a, b and c are determined by a certain 


linear system of equations, (see [ Ref. 6]). 


C. NIELSON“FRANKE QUADRATIC TRIANGLE METHOD 

All methods listed under triangle based blending methods 
[ Ref. a] are conceptually similar to the inverse distance 
weighted methods. The interpolation function for the quad- 
ratic triangle method is 


P (,Z ) = Ws (X,Y) "Q, (X,Y) i Ws (X,Y) "Qs , ) rs Wy, (X,Y) “Oy (X,Y) 


The first step in the present method is to partition the plane 
into triangles by connecting neighboring data points based 
upon the min-max angle criterion as described by Franke in 

[ Ref. 5 |. The quadratic triangle method uses the inverse 
distance weighted quadratic Q, (X,Y) which is also used in the 


modified quadratic Shepard's method, see section 2.A. 
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A significant difference from the inverse distance weighted 
method is that the weight functions for the triangle based 
blending method are based on a triangulation of the convex 
hull of the point set { (x, +%) - The method is not sensi- 
tive to the triangulation technique. If a triangulation is in 
existence for other purposes, it can be used. 

In short, the method is to define the nodal functions 
Q.rk=1,..-5N, as in the modified quadratic Shepard's method, 
form a triangulation of the points, determine the weight 


functions [ Ref. = | and eGempute interpolant F(X,Y). 


D. NIELSON'S MINIMUM NORM NETWORK METHOD 

Nielson's minimum norm network method | Ref. 15], can 
be found under the finite element based methods in Table l. 
Like the Nielson-Franke quadratic triangle method, this global 
method is based upon a triangulation. 

The method consists of three separate steps; triangulation 
curve network, and blending | Ref. 7]. The points 
),kK=1,...,N are used as the vertices of a trlangulation 


Kk 
of the convex hull of these points. Alternatively, if a tri- 


(XY 


aggulation already exists, it can be used. The curve network 
step involves the solution of a certain minimum norm problem 
and requires first partial derivatives in its discretized 

form. These are obtained by assuming a cubic variation along 
each edge in the triangulation and minimizing the integral of 


the second derivative squared [ Ref. 7]. 


ie 





This method does not provide extrapolation outside the con- 
vex hull, but the triangular blending method can be used to 
em@eenad the curve network to the entire hull. In a triangular 
exterior region, the function is taken to be the linear 
Function determined by the value and slopes at the vertex 


[ Ref. 7 |. 


E. GENERALIZED NEWTON DELTA SUM BICUBIC SPLINE METHOD 

Foley's methods [ Ref. 4 | involve several ideas. The 
interpolant is taken as either the generalized Newton inter- 
polant, or a form of Shepard's method. The use of generalized 
Newton type interpolant is involved in them prominently, 
because the best performance for various test functions is 
provided by the interated delta sum methods using the generalized 
Newton interpolant with natural bicubic splines. 


The generalized Newton interpolant takes the form 





N 
[ee Ces 5 al W (X,Y) , 
N k=1 * 
EE f= Tt (aay) 
ie k-l &k ae 
where a= 
and W (X,Y) has the property eee) 100) es ee ae oe 


kk 
This function is dependent upon the order of the points, and 


So Foley's scheme involves an ordering process [ Ref. 4]. 
After constructing a bicubic spline interpolant for the grid 
points, one adds a generalized Newton interpolant for the 


difference between data and the bicubic splines, obtaining 


20 





an interpolant. This process is termed a 'delta sum' by 


Foley [ref. 4]. 


Pe. BARDY'S MULTIQUADRIC METHOD 
Hardy's multiquadric method [ Ref. 8,9 | can be thought 
of as being under the global basis function type. A multi- 
quadric surface can be represented by 
IY? 


N 2 2 
i C. (Xi K.) + (Y -Y¥Y ) +€¢C =—oe  -L=ljcselsN % 
3) Dee 


in 
The given problem data provide a set of cartesian coordinates 
on the surface ranging from XY 72, to Koel oui the quad- 
ratic term coefficients C, to Cy are unknown. The coordinates 
of the N data points can be substituted into the definition 
of the surface and only the coefficients C, are left as 

J 
unknown. A problem here is that of solving a system of N 
linear equations with N unknowns. A nominal value for C was 
determined by computational experiment and the value C = 0 


was chosen. 


ST 

“2 
If a vector of unknowns is to be [c, ] = : | 
used, cGefine ; 

On 


then each known element becomes 


Ze 





Zz Z 
OS + (Y -¥,) =o, 
ee 7 ess 


and the matrix [ a, . | = A, which is (NXN) coefficient matrix. 


Upon 
ra] 
. 1| 
Zo) 
defining which pz | = . 
4 
se 2) 
A multiguadric surface reduces to AxC = Z and, as usual, the 


solution is C = avtg. 


When the Known constraints . are substituted into the 
multiquadric surface formula, we have the required equation 
of a surface, which fits the data points exactly and which 


provides logical interpolation at intermediate points. 
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IIIT. COMPARISON OF INTERPOLATION METHODS USING 
RALSTON VALLEY DATA 

Contour construction for seismic p-wave velocities and 
depth data can be seen in Appendix C and in Appendix F, 
respectively, for the various interpolation methods and with 
the three different contour plot techniques. Each method 
appearsin the order given in Table 1. Each method is dis- 
played uSing first the subroutine CONTUR, then the subroutine 
CONISD, and lastly the subroutine PLT3D1l. Subroutine CONTUR 
Gaamot outline the Ralston Valley, but gives contour con- 
Structions in rectangular form with the x-axis ranging from 
one to seventeen miles and the y-axis from one to twenty five 
Miles. The subroutine CONISD can outline the Ralston Valley 
with boundary lines. All three dimensional pictures were 
Generated by subroutine PLT3Dl. 

Appendix D shows the mean and standard deviation of the 
various contour constructions for the seismic p-wave velo- 
cities. Numerical values for the mean and standard deviation 
at each grid point come from six different interpolation 
methods. Appendix G is similar to Appendix D, but it is for 
the depth of the surface soil layer. Convex contour lines 
are formed for the standard deviations values as can be seen 
in Appendices D and G. Small standard deviation values 


appear in data point region, and become larger as one moves 


Z5 





toward the edges. In other words, more accurate prediction 
occurs in data point region and less accurate prediction out- 
sede this region. 

Deviations from mean values for all interpolation methods 
can be seen in Appendix E for seismic p-wave velocities and 
in Appendix H for depth values. The 1, 10, 50 and 100 contours 
show the absolute differences between the given method and 
mean values for seismic p-wave velocities. The 0.1, 0.5, 1.0 
contours for depth data display quite convex regions for all 
methods. Clearly, it can be seen that; there are no signifi- 
cant differences between methods in terms of deviations from 
mean values for grid points especially in data region. This 
result shows that in the data points region all methods work 
with almost same accuracy and outside the region there are 
too great differences. There is no way to compare them out- 


Side the region. 
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ive CONROUR PLOT  TECHNLOUES 


The interpolation values provided by any of the methods 
can be seen in the computer outputs, but their interpretation 
mo egnly difficult without contour construction. For inter- 
pretation purposes, computer programs for drawing plots on 
an off-line plotter are available at the Naval Postgraduate 
School computer center. 

The subprograms CONTUR, CONISD and PLT3Dl serve the com- 
mon purpose but have different capabilities. 

All the subprograms have the same limitation in that they 
use a linear interpolation process to find the values of the 
points along contour segments. Also, there is a limit in the 
resolution of the off-line plotter, namely, .005 inch in both 


the X and Y directions for the versatec plotter. 


& SUBROUTINE CONTUR 

Given a matrix of numerical values of Z = (X,Y), the 
program CONTUR generates a contour graph on which one or 
more contour levels are drawn. This non-IMSL subroutine 
uses a two-dimensional array and turns out a two-dimensional 
off-line plot. Interior and exterior contour segments and 
labels (if requested) for exterior contour segments can be 
seen. Labeling of the interior segments represent local 


maxima. 
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1) 


2) 


1) 


2) 


4) 


>) 


6) 


Advantages: 

The subroutine has a minimum number of arguments (10) 
compared with other plot techniques. 

A rather big two-dimensional array can be used, up to 
ms? data pomnts. 

It takes approximately one fiftieth of the CPU time that 
would be required by subroutine CONISD. 

Disadvantages: 

Irregularly spaced data cannot be contoured. 

Irregular boundary capability is not available, and there 
1S no option for a 'cut-out' area. 

The width and the height of the contour graph must be 
specified as integers, and are not adjustable scale 
parameters. 

The scale values (1f requested) one inch apart on the 
exterior frame of the contour graph increase from north 
to south. This 1s not oriented properly for our use, 
which is a lower left grid oriented system. 

Limited labeling. The x-axis and y-axis cannot be 
labeled. 


There is no option for smoothing the contour lines. 


SUBROUTINE CONISD 


This non-IMSL subroutine produces a drawn contour map on 


an off-line plotter for given irregularly spaced data points. 


Each data point is a triad of X,Y and Z values where 


2 = 


(X,Y). There may be one or more ‘cut-out areas' which 


ame not COMNTOUr. 
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Advantages: 

1) tIrregularly spaced data can be drawn. 

a trzegular boundary capability is available for one or 
more ‘cut-out areas’. 

3) The scale for the x-axis and y-axis are defined as real 
numbers, 1.e., adjustable scale. 

4) X-axis and y-axis can be labeled. 

5) There is an option for smoothing the contour lines. 
Disadvantages: 


1) It takes much more CPU time than the other plot techniques. 


e. SUBROUTINE PLT3D1 

A subroutine PLT3Dl produces three-dimensional perspec- 
tive or isometric plots on an off-line plotter for given a 
two-dimensional matrix of numerical values of Z = (X,Y). 
Labeling and scaling are not available for this reason, it 
can be used only for general purposes. 

A memorandum about this subprogram [ Ref. 12 | can be 
requested from W.R. Church computer center. This memorandum 
includes a listing of the computer program, the algorithm 
that is used by the program, and the axes orientation, rota- 
tion and projection are clarified by illustrations. 

Advantages: 

1) This routine gives a nice overview to given area by 

aS Sey) . 

2) The vantage point can be changed. 


3) The scale is adjustable. 


ae 





Disadvantages: 
i terme 1S mo option for ‘cut-out areas’. 


2) There is no labeling. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


Selsmic p-wave velocities or depth values for any grid 
point in Ralston Valley can be found from computer outputs or 
contour constructions for any particular interpolation method 
Or for the mean of the six different methods. The standard 
deviations for each grid point are also available. The means 
and standard deviations may be useful for the comparison of 
methods. 

The computer outputs for standard deviations in Appendix 
E and H show that accuracy 1s high in the data point region 
and at some undecided distance from this region; but this 
distance cannot be measured. In the data point region all 
interpolation methods gives essentially the same result, 
regardless of whether they are global or local. 

Outside the data point region all methods generate rather 
different values and there is no way to decide which one is 
close to true values. If the point of interest is in the data 
point region any of the interpolation methods can be used. 

According to Appendix C and F; none of the methods can 
find a maximum or minimum point outside the region. Table 2 
gives approximate CPU times for all of the studied inter- 
polation methods. If minimal computer time is required for a 


given problem, the triangle based quadratic method is suitable. 
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